---
title: "Class 6"
author: "Colin Kuehl"
date: "`r format(Sys.time(), '%B %d, %Y')`"
output: pdf_document
editor_options: 
  chunk_output_type: console
---

```{r,echo=FALSE,message=FALSE,warning=FALSE}
library(formatR)#may need to install this first
# Set so that long lines in R will be wrapped:
knitr::opts_chunk$set(tidy.opts=list(width.cutoff=80), tidy=TRUE)
```

```{r}
setwd("~/Dropbox/POLS 641/CMU Su26/cmucourse26/ClassCode/Day 6")

```

# Many Means

We'll start out illustrating concepts from discussion using R

## Distribution of the Mean

```{r}
getmean=function(N, mean, sd){
  opinion=rnorm(n = N, mean=mean,sd=sd)
  meanopinion=mean(opinion)
  return(meanopinion)} 
#Building a function in R. This is beyond the scope of the class

manymeansa=replicate(n = 10, getmean(20, mean=65, sd=6))
hist(manymeansa, main="Many Means 10 Draws", sub=mean(manymeansa))
mean(manymeansa)
abline(v=mean(manymeansa), col="orange")

#look with different n
manymeansb=replicate(n = 20, getmean(20, mean=65, sd=6))
hist(manymeansb, main="Many Means 20 Draws", sub=mean(manymeansb))
mean(manymeansb)
abline(v=mean(manymeansb), col="orange")


manymeansc=replicate(n = 100, getmean(20, mean=65, sd=6))
hist(manymeansc, main="Many Means 100 Draws", sub=mean(manymeansc))
mean(manymeansc)
abline(v=mean(manymeansc), col="orange")


manymeansd=replicate(n = 1000, getmean(20, mean=65, sd=6))
hist(manymeansd, main="Many Means 1k Draws", sub=mean(manymeansd))
mean(manymeansd)
abline(v=mean(manymeansd), col="orange")

manymeanse=replicate(n = 1000, getmean(20, mean=65, sd=6))
hist(manymeanse, main="Many Means 1000 Draws", sub=mean(manymeanse))
mean(manymeanse)
abline(v=mean(manymeanse), col="orange")

manymeansf=replicate(n = 10000, getmean(20, mean=65, sd=6))
hist(manymeansf, main="Many Means 10k Draws", sub=mean(manymeansf))
mean(manymeansf)
abline(v=mean(manymeansf), col="orange")

#Larger sample sizes should be closer to "true mean" but still variation
abline(v=65, col="blue")
abline(v=mean(manymeansa), col="darkorange4")
abline(v=mean(manymeansb), col="darkorange3")
abline(v=mean(manymeansc), col="darkorange2")
abline(v=mean(manymeansd), col="darkorange1")
abline(v=mean(manymeanse), col="darkorange")
abline(v=65, col="blue")



mean(manymeansa)
mean(manymeansb)
mean(manymeansc)
mean(manymeansd)
mean(manymeanse)
mean(manymeansf) #as n increases will get closer to true mean

sd(manymeansa) # . . . and sd
sd(manymeansb)
sd(manymeansc)
sd(manymeansd)
sd(manymeanse)
sd(manymeansf)
```

################################################################################################# 

#Bechdel Test

Load the data: These data comes from this [fivethirtyeight article](https://fivethirtyeight.com/features/the-dollar-and-cents-case-against-hollywoods-exclusion-of-women/)[archived here](https://web.archive.org/web/20260103173002/https://fivethirtyeight.com/features/the-dollar-and-cents-case-against-hollywoods-exclusion-of-women/)

```{r}
load("bech.RData")
```

Load Packages: we'll use ggplot briefly again today(we dont' need to install package again)

```{r}
library(ggplot2)
```

Find the number of variables and observations - *What is our unit of analysis?*

```{R}
head(bech)
```

VARIABLE NAMES

```{r}
names(bech) 
variable.names(bech) 
```

What time range does the data cover?

```{r}

```

## Exploring the data

```{r}
table(bech$clean_test)

plot(bech$domgross, bech$intgross)

prop.table(table(bech$binary))

prop.table(table(bech$clean_test, bech$binary))

barplot(table(bech$clean_test)) #show ugly bar plot

```

Bechdel said that if a movie can satisfy three criteria it passes the test 1. there are at least two named women in the picture 2. they have a conversation with each other at some point 3. that conversation isn’t about a male character

FIND THE RANGE OF GROSS REVENUE(we'll just look at domestic and standardize to 2013 dollar value)

```{r}
range(bech$domgross_2013)
max(bech$domgross_2013) #Bonus: What movie is this? Hint: View(yourdata)
min(bech$domgross_2013) #Cheers to you if you've actually watched this
```

FIND STANDARD DEVIATION AND VARIANCE - Formula for Standard Deviation: $$\sigma = \sqrt{\frac{\sum\limits_{i=1}^{n} \left(x_{i} - \bar{x}\right)^{2}} {n-1}}$$ More info on using math symbols in Rmarkdown is here: <https://rpruim.github.io/s341/S19/from-class/MathinRmd.html>

```{r}
sqrt(sum((bech$domgross_2013-mean(bech$domgross_2013))^2/(length(bech$domgross_2013)-1))) #How would we turn this into the variance?
```

Of course, it is much easier to just have R do it

```{r}
sd(bech$domgross_2013)
var(bech$domgross_2013)
sd(bech$domgross_2013)^2 #Variance = standard deviation squared.
sqrt(var(bech$domgross_2013)) #Standard deviation = square root of variance.
```

TAKE A LOOK AT THE DISTRIBUTION OF DOMESTIC GROSS

```{r}
hist(bech$domgross_2013, breaks = 15) #breaks=  number of bars
hist(bech$domgross_2013, main="Distribution of Gross Revenue", xlab="Gross Revenue", col="blue",border="green",breaks=30, las=2) #bedazzled
```

LOOK AT THE DENSITY PLOT OF DOMESTIC GROSS

```{r}
plot(density(bech$domgross_2013),main="Density Distribution of Gross Revenue")

ggplot(bech, aes(x=domgross_2013))+geom_density()+
  labs(x="Domestic Gross 2013", y="Density", title="Density Distribution of Gross Revenue", caption="ggplot version")
```

How are histogram and density plot different? Which is closer to a pdf?

```{r}
hist(bech$domgross_2013, main="Histogram of of Gross Revenue",
xlab="Domestic Gross 2013",breaks=20)

ggplot(bech, aes(x=domgross_2013))+geom_histogram()+
  labs(x="Domestic Gross 2013", y="Density", title="Histogram of Gross Revenue", caption="ggplot version")
```

## MAKE A BOXPLOT

```{r}
boxplot(bech$domgross_2013, main="Boxplot of Gross Revenue", col="cornflowerblue")
abline(h=mean(bech$domgross_2013), col="deeppink", lwd=2)
```

A boxplot shows the median, 25% and 75% quartile, and important outliers

ADD LINES TO INDICATE THE MEAN AND MEDIAN, AND MAKE A LEGEND *(Make sure all the code stays together)*

```{r}
plot(density(bech$domgross_2013), main="Density Plot of Gross Revenue", xlab="Gross")
abline(v=mean(bech$domgross_2013), lty=2, lwd=2, col="red") # put a vertical line at the mean 
abline(v=median(bech$domgross_2013), lty=4, lwd=2, col="blue")
legend("topright", legend=c("mean", "median"), lty=c(2,4), col=c("red","blue"))

options(scipen=10)#This stops use of scientific notation
```

## Take Random Samples and see how mean and median changes

```{r}
sample1 <- bech[sample(nrow(bech), 100), ] #Randomly samples 100 movies

plot(density(sample1$domgross_2013), main="Density Plot of Gross Revenue", xlab="Gross")
abline(v=mean(sample1$domgross_2013), lty=2, lwd=2, col="red") 
abline(v=median(sample1$domgross_2013), lty=4, lwd=2, col="blue")
legend("topright", legend=c(mean(sample1$domgross_2013), median(sample1$domgross_2013)), lty=c(2,4), col=c("red","blue"))


sample2 <- bech[sample(nrow(bech), 100), ]
plot(density(sample2$domgross_2013), main="Density Plot of Gross Revenue", xlab="Gross")
abline(v=mean(sample2$domgross_2013), lty=2, lwd=2, col="red") 
abline(v=median(sample2$domgross_2013), lty=4, lwd=2, col="blue")
legend("topright", legend=c(mean(sample2$domgross_2013), median(sample2$domgross_2013)), lty=c(2,4), col=c("red","blue"))


mean(sample1$domgross_2013)
mean(sample2$domgross_2013)

median(sample1$domgross_2013)
median(sample2$domgross_2013)
```

LOOK AT THE SAMPLING DISTRIBUTION OF REVENUE AND BUDGET *(Let's put the four figures arranged in 2 rows and 2 columns)*

```{r}
par(mfrow=c(2,2)) ###This line of code tells R you want a 2 x 2 of smaller graphs
hist(bech$domgross_2013, main="Histogram of Domestic Gross", breaks = 15)
plot(density(bech$domgross_2013), main="Distribution of Domestic Gross")
hist(bech$budget_2013, main="Histogram of Budget", breaks = 15)
plot(density(bech$budget_2013), main="Distribution of Budget")
```

THIS IS NOT A NORMAL DISTRIBUTION SO MAY WANT TO CREATE LOGGED VERSION We will talk about transforming variables later

```{r}
par(mfrow=c(1,1)) #Return to just looking at one graph at a time
bech$logdomgross13 <- log(bech$domgross_2013)#creates new logged variable
hist(bech$logdomgross13, main="Logged Domestic Revenue")#Looks more normal, harder to interpret
boxplot(bech$logdomgross13, main="Logged Domestic Revenue")

hist(bech$domgross_2013, main="Domestic Revenue")#Looks more normal, harder to interpret
boxplot(bech$logdomgross13, main="Logged Domestic Revenue")

```

## LET'S USE STATS TO MAKE COMPARISONS

First, find averages for films that pass and fail the Bechdel test (there are several ways to do this).

```{r}
pass <- bech[which(bech$binary=="PASS"), ] #Creates a new dataset of only those that pass the bechdel test
fail <- bech[which(bech$binary=="FAIL"), ]
mean(pass$domgross_2013)
mean(fail$domgross_2013)
#How many films passed and how many failed?

mean(pass$budget_2013)
mean(fail$budget_2013)

ggplot(bech, aes(x=logdomgross13, fill=binary)) +
    geom_density( alpha=.5, position="identity") #overlay both distributions
```

Second, simple algebra to find difference.

```{r}
grossdiff <- mean(pass$domgross_2013)-mean(fail$domgross_2013)
grossdiff
```

*Another method for finding mean of subset (more here: <http://www.statmethods.net/management/subset.html>)*

```{r}
xbar_pass_t <- mean(bech$domgross_2013[bech$binary=="PASS"])
xbar_fail_t <- mean(bech$domgross_2013[bech$binary=="FAIL"])
DIFF <- (xbar_pass_t - xbar_fail_t)
DIFF
```

Show Both Box Plots

```{r}
boxplot(pass$domgross_2013)
boxplot(fail$domgross_2013)
#or from original dataset
boxplot(bech$domgross_2013~bech$binary) 
#using ggplot
ggplot(bech, aes(binary,domgross_2013))+geom_boxplot()+
  labs(y="Domestic Gross 2013", x="Bechdel Test") #creates boxplots based on values of binary

#to illustrate we can use this for variety of categories(and logged)
ggplot(bech, aes(test,logdomgross13))+geom_boxplot()+
  labs(y="Domestic Gross 2013", x="Bechdel Test")+ theme(axis.text.x=element_text(angle = -45, hjust = 0))



ggplot(bech, aes(as.factor(year),logdomgross13))+geom_boxplot()+
  labs(y="Domestic Gross 2013", x="Bechdel Test")+ theme(axis.text.x=element_text(angle = -90, hjust = 0))

```

## FINDING RELATIONSHIPS *(Let's see if bigger budgets are related to a higher gross)*

```{r}
par(mfrow=c(1,1)) #show one graph at a time
plot(bech$budget, bech$domgross_2013, main="Gross Revenue and Budget", xlab = "Budget", ylab = "Gross Revenue")
```


ggplot version
```{r}
ggplot(bech, aes(x = budget, y = domgross_2013)) +
  geom_point() +
  labs(title = "Gross Revenue and Budget",
       x = "Budget",
       y = "Gross Revenue") # simple version

ggplot(bech, aes(x = budget, y = domgross_2013, color = binary)) +
  geom_point(alpha = 0.6) +
  scale_color_manual(values = c("FAIL" = "firebrick", "PASS" = "steelblue")) +
  labs(title = "Gross Revenue and Budget",
       x = "Budget", y = "Gross Revenue", color = "Bechdel Test") +
  theme_minimal()
```




IT'S HARD TO SEE A PATTERN FROM THE PLOT, SO LET'S MAKE A LINEAR REGRESSION MODEL: REMEMBER THE MODEL IS lm(Y \~ X, data=nameofdata)

```{r}
plot(bech$budget, bech$domgross_2013, main="Gross Revenue and Budget", xlab = "Budget", ylab = "Gross Revenue")
model1 <- lm(domgross_2013~budget, data=bech) #Making an object for a linear regression model. Notice how I access the variables used... 
summary(model1) #regression output. We'll spend lots of time here tomorrow
abline(model1, col = "purple", lty=2) #This puts the trend line

abline(lm(pass$domgross_2013~pass$budget_2013), lty=3, lwd=2, col="pink") #for fun add lines for films that pass bechdel test
abline(lm(fail$domgross_2013~fail$budget_2013), lty=4, lwd=2, col="blue")
legend("topright", legend=c("All", "Pass", "Fail"), lty=c(2,3,4), col=c("purple","pink", "blue"))
```

What do the different slopes mean?What does this tell us about the profitability of inclusing in the movie industry?

Another Subset: We may also want to exclude outliers

```{r}
middle <- subset(bech, budget_2013 >= 1000000 & budget_2013 < 10000000) #only budgets between 1 and 10 millsion
summary(middle$budget_2013)

#do the relationships change when we control for outliers(hint: use copy and paste plus find and replace)



#What about when we use the log of gross



```

################################################################ 
